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Abstract. This work proposes a parametric equilibrium distribution function !F eq to 
be applied to the gyrokinetic studies of the Finite Orbit Width behavior of guiding 
centers representing several species encountered in axisymmetric tokamak plasma, as 
fusion products, thermal bulk and energetic particles from Ion Cyclotron Radiation 
Heating and Negative Neutral Beam Injections. 

After the analysis of the basic results of orbit theory obtained with a particularly 
convenient orbit coordinates set, it is shown how the proposed T eq satisfies the two 
conditions that make it an equilibrium distribution function: (i) it must depend only 
on the constants of motion and adiabatic invariants, and (ii) the guiding centers must 
remain confined for suitably long time. 

Furthermore, the T eq can be modeled, with a proper choice of its parameters, to 
reproduce the most common distribution functions. A local Maxwellian distribution 
function is obtained for the thermal plasma in the Zero Orbit Width approximation. 
For the fusion a particles, F eq can also reproduce the Slowing Down (SD) distribution 
function. More generally, for supra-thermal particles, when external heatings are 
present, such as (N)NBI and ICRH, the proposed model distribution function shows 
similarities with the anisotropic SD and the biMaxwellian distribution functions. 

T eq can be used to fit experimental profiles and it could provide a useful tool for 
experimental and numerical data analysis. Moreover, it could help to develop analytical 
computations for facilitating data interpretation in the light of theoretical models. This 
distribution function can be easily implemented in gyrokinetic codes, where it can be 
used to simulate plasma also in the presence of external heating sources. 

PACS numbers: 52.30.Gz, 52.20.Dq, 52.25.-b, 52.50.-b, 52.20.-j 
Keywords: Equilibrium Distribution Function, Orbit Theory, Plasma Heatings, 
Gyrokinetics. 

1. Introduction 



Gyrokinetic theory in general and gyrokinetic simulations in particular make often use of 
an initial distribution function of guiding centers (GCs), usually indicated by F . Initial 
distribution function must represent a slowly evolving equilibrium for a perturbative 
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approach. There are two conditions that F must satisfy in order to be considered an 
equilibrium T eq one: 

it must depend only on the invariants of motion (1) 

and 

the GCs must remain confined for suitably long time. (2) 

The meaning of the first sentence is clear: the function of constants is constant. The 
second condition takes into account that laboratory plasma is placed into a limited 
portion of space; if the particles are contained within this region only for a finite time 
interval and then lost, consequently the corresponding particles distribution will be 
impoverished in the protracted time. 

The work is set up to progressively build the equilibrium distribution function. In 
the next section, it is shown the chosen set of constants of motion (COM) in order to 
describe the unperturbed orbits characterizing the equilibrium. There is no need of 
entering into the details of the Guiding Center (GC) transformation. The reader who 
is not confident with this coordinate transformation can consult [1] or one of the many 
introductory articles on this topic. Neverthless, it is explained the physical assumptions 
needed to consider as constants the following GC quantities: (corresponding to 
the axisymmetry angular momentum), the kinetic energy per unit mass w and the 
generalized pitch angle variable A, defined as the ratio between the magnetic moment \x 
and w. In the text, the set (V<f,,w, A) will be referred to as the Quasi Invariants (QIs) 
set. The Section (2.2) provides a brief outline on how it is currently solved the problem 
of assigning an equilibrium distribution function in gyrokinetic codes [2, 3, 4, 5, 6, 7, 8]. 
The difficulties arising from that common procedure will be the starting point for the 
developing of an alternative approach. Section 3 addresses the above condition (2), 
where a GC is considered lost when a portion of its unperturbed orbit will be out of the 
plasma section. The main topic of this Section is devoted to specify the topology of the 
projected orbits onto the poloidal section. Indeed, to better understand an equilibrium 
distribution function of QIs, it is convenient to consider it as a distribution of (projected) 
orbits. It will be useful to know if GCs are, not only confined or lost, but also, trapped 
or passing. 

In this analysis, the poloidal magnetic flux coordinate ip, in addition to V^, w and 
A, are used as orbit coordinates. In the light of this set of coordinates, it will be 
convenient to show the basic results of the orbit theory [9, 10, 11, 12, 13, 14, 15], such 
as the orbits classification, the derivation of a graphical method to promptly deduce 
the orbit shape, the definition of the orbit average, the expression of the characteristic 
frequencies (transit and bounce frequencies) and the expression of the second invariant. 

Finally, in Section 4 the parametric T eq equilibrium distribution function is built 
adopting the following guidelines: it must have Boltzmann-like behavior, the energy 
must be represented in terms of QIs and it must be mathematically tractable. Here, the 
similarities found in various limiting behaviors are represented: when the proposed J-" eg 
is compared with the commonly used distribution functions in tokamak plasma physic 
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such as the local Maxwellian, the Slowing Down (SD), the biMaxwellian distribution 
function and so on. 

Section 5 qualitatively describes the various terms that constitute the obtained T eq . 

An appendix has been inserted solely to support the reader in order to easily refer 
to some known formulas for the description of the magnetic field with the adopted 
convention. Appendix A briefly considers two particularly used coordinates reference 
systems: magnetic flux coordinates and Shafranov coordinates. The Shafranov geometry, 
characterized by circular plasma poloidal section, is the example geometry adopted to 
visualize some results useful for more general plasma poloidal sections. 

In this work the vectors are not indicated with bold letters. As for example, the 
guiding center velocity V, the electric field E, the magnetic field B, the unit vector b in 
the direction of the magnetic field, the perpendicular (to the magnetic field) velocity v±, 
the drift velocity vd, and the grad operator V are vector quantities, whilst the parallel 
(to the magnetic field) velocity v\\ G [—00, +00], the velocity magnitude v and the radial 
coordinate r are scalars. The symmetry axis component of the angular momentum 
Lz and the toroidal component of the canonical angular momentum V$ are vector 
components. When there is a possibility of confusion it will be specified the scalar or 
vector character of the adopted symbol. 

Moreover, natural units (n.u.) are employed: the speed of light c = 1. 

2. Preliminary concepts 

It is well known that the particle energy £ = m s v 2 /2 + q s A (where A is the electric 
potential and m s ,q s specify the considered particle species with its mass and its 
charge) and the (canonical) angular momentum component along the symmetry axis 
Lz = m s Rv<f, + qsRAtf, (where and are respectively the toroidal component of 
the magnetic vector potential A and of the particle velocity v) are COMs, assuming 
a charged particle (m s , q s ) in non relativistic regime, the toroidal symmetry and the 
presence of the only Electro-Magnetic field. Lz is frequently named canonical toroidal 
momentum. Finally the magnetic moment p = v"j_/(2\B\) is an adiabatic invariant. 

These particle invariants undergo the Guiding Center (GC) transformationf which 
reduces the dimensionality of the velocity space from 3D to 2D, removing the gyroangle 
7, the angle that fixes the direction of v±, up to a wanted order in 8 = p/L (where p is 
the Larmor radius and L is a typical scale length of the system). As for the gyroangle, 
also the conservation property of the transformed GC COM is usually verified up to a 
given order in 5. Thus, preserving the same symbols and fixing a given 8 order, Lz,£ 
and p are (approximate) GC COMs. 

Here a slightly different viewpoint is adopted: the constants of motion in GC 
coordinates are exact constants, whilst the corresponding quantities, expressed in 
particle coordinates, become approximate. In a schematic manner, it is common practice 

f The reason why it is used the GC instead of the gyrocenter transformation is because it is analyzed 
the equilibrium, where the fields are unperturbed. 
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considering 

(exact) COM in particle coord. = (approximate) COM in GC coord. + 0(S n ) 
while in this paper 

(exact) COM in GC coord. = (approximate) COM in particle coord. + 0(5 n ). 

The results must be correct up to the (n-l)-th 5 order for both descriptions. Although 
the first and most used approach is more realistic, the second one is preferable when 
there is no need to come back to the particle coordinates once arrived at the GC 
description. Ref.[16] uses the same approach referring to COM in unperturbed GC 
dynamical reduction system. The use of the second approach is clearer than the first 
one because it fulfills the exact conservation of some GC quantities. The imposed 
constants of motion for the GC description are renamed QI to differentiate the two 
approaches: Quasi Invariants are exact Constants Of Motions in GC coordinates. 



2.1. Quasi Invariants 

The unperturbed GC projected orbit onto a poloidal section will be described in the 
toroidal coordinates r, 9 and 0, by the following relations: 

f — V • Vr, 

3 

9 = V ■ V0, 

while the toroidal motion is described by: 

= V ■ V0, (4) 

where V is the GC velocity. Here, the definition of the poloidal angle 9 is quite arbitrary: 
it is a 27r periodic coordinate and it lies in a plane orthogonal to the toroidal unit vector 
e<£. Even the definition of the radial coordinate can differ but it is required labeling the 
magnetic flux: r = r(ip). 

When the E x B drift, the V-B drift and the curvature drift is retained correctly 
whilst the 0(<5 2 ) drift is neglected, then V can be expressed as: 



where b 



V = Qv\\b + v D = Q 
B/\B\ 



m s v\\ „ 

V\\b+ ~4^V X 



>ll&) 



(5) 



Qs\B 

Vd is the drift velocity and Q — 1 + 0(5) depends on how 
have been defined v\\ in the GC transformation: commonly v\\ = V ■ b and Q = 
1 + m s /(q s \B\)b ■ V x (v\\b) as in [17, 18], otherwise if v\\ = (V - v D ) ■ b then Q = 1 
recovering the original expression proposed by [19]. 
The equations in (3) are explicited as follows: 



f = Q^JrV x (u,|6) • Vr 
Qs\B\ 

Tfi V\\ 

9 = Q Vll b-V9 + — ^[V x (u,|6) • V9 
Qs\B\ 



(6) 
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For convenience, the magnetic field is written as follows (A. 13): 

B = Vip x V0 + F(ip) V0, (7) 

where the radial component of the plasma current density is required to be zero to 
ensure force balance (otherwise F(ip) should be replaced by F{ip,9), see (A. 11)). 

The terms b ■ V0, V x (v\\b) ■ V6> and V x (v\\b) ■ Vr, thanks to the symmetry, 
becomes: 

V x („„&) • W = V(^) xVf V^ = - -Lfl^(^) (9) 

V x (u||6) • Vr = V(^) xVf Vr= ^-d e {^-) (10) 

l"l \9 \B\ 

where y/g = (Vr x V6 • V0)" 1 is the spatial Jacobian for the toroidal coordinates 
transformation, ip = ip(r) is invertiblej and the prime indicates the radial derivative 
with ip' < 0. Because V9 together with Vr is orthogonal to V0 = e^/R, then 
V x (r>|j \B\~ l Vi> x V0) is parallel to (it does not give any contribution if scalarly 
multiplied per W or Vr). From (6,8,9,10) and with f = ip/ip' and F = m s F/q s , the 
following relations are obtained: 

= ^S-^L^i = - J. (u) 

I-BIV9 l B l </>' + 8 r (F«n/|B|) 1 

The constancy of V$ defined as V<$> = ip + Ff||/|_B| is shown putting ip on the RHS, 
changing the sign and multiplying for the denominator: 

o = * + ^ + i^ = 4 + ±$ = ±« + $ ) = i>, («) 

The statement in the opposite direction is not true: = does not imply the 
drift velocity in (5). Indeed, the constancy of can also take into account a toroidal 
flow: V = Q v\\b+ (m s v\\) / (q s \B\)V x (v\\b) +lZ(r,0)V(j), where 1Z stands for rotation. 
This is because the equation (4) has not been considered yet. What has been shown is 
a clear correspondence between V<p and the expression of the drift velocity (5) due to 
the toroidal symmetry. 

When if) = RA^ (A. 9) is substituted within V^, it becomes clear also the reason 
why the gyrokinetic community refers to V<p as the canonical toroidal momentum: 

+ (13) 
uj c q s 

for uj c the cyclotron frequency. 

% For simplicity, the magnetic flux surfaces are considered nested. 
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The GC motion is further simplified by requiring the constancy of w and ji. It 
is worth noting that the condition on w is coherent with the neglecting of the electric 
potential. Indeed in drift ordering Aq behaves as 0(5): 

w=^ + fi\B\^£; (14) 

however, the present analysis could straightforward include an electric potential if it 
depends solely on ip. 

From w = and fi — 0, the time derivative of v\\ becomes [16]: 

m = -V-V-V\B\. (15) 

v \\ 

Now there is a clear correspondence between (3) and (15) with the following system of 
equations: 

w = 0, 

A =0, (16) 
^ = 0. 

The GC projected orbit may be defined by the initial conditions: w = wo, A = Ao and 
Vcj, = "P</,o§- As mentioned above, QIs are the imposed constants w, A and V^. It becomes 
clear that a distribution function depending on the QIs will describe a distribution of 
GC (projected) orbits. 

When the orbits behavior is considered in a simulation or in a theoretical analysis, 
it is necessary to take into account the finite orbit width (FOW) effects. It is possible to 
estimate the relevance of the FOW directly from the expression of V^. Indeed, the term 
Fv\\/oo c in (13) expresses how far apart an orbit will be from the poloidal flux surface 
coordinate ip = e.g. characteristic banana orbits describing trapped particles will 
have the tip of the banana at ip — because here v\\ = 0, whilst the banana orbit 
width will depend on the maximum and minimum reachable values of v\\. 

Several possibilities may arise. Concerning the electrons, u c is big enough. In this 
case it is possible to describe the electron orbits directly on the flux surface, because 
■0 ~ This is called the small orbit width (SOW) case. If oj c — > oo then if) = V^. This 
is called the ideal ZOW (Zero Orbit Width) case and the equilibrium can be described 
by a distribution function depending on (ip,w,\). In the ZOW case, the often used 
equilibrium distribution function is the local Maxwellian distribution function (where e 
is for electrons): 

f M U w ) = e -»HW (17) 



where v te = \JT e /m e . A useful extension of this distribution function, taking into 
account SOW effects, is depicted by: 

jM{r<f>,W) (2^)3/2^3^) > ^) 

§ Sometimes the sign of vu has to be specified, how it will be clarified in Section 3. 
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known as canonical Maxwellian distribution function [2]. The local Maxwellian 
distribution function can be considered the ZOW limit of the canonical Maxwellian 
distribution function, appropriate for SOW effects. In the following some difficulties, 
arising specially when large orbit width effects are not negligible, will be described. 

2.2. The current way to describe a gyrokinetic equilibrium with Finite Orbit Width 
effects 

Several difficulties arise when the considered orbit has a large width. Firstly it is 
analyzed what happens when the approximation if) ~ V<j, fails. This is more evident 
for the passing particles, for which v\\ ^ almost always, causing a shift to the 
whole orbit respect to the flux surface value ip = V^, as can be easily visualized 
anticipating the Figure 1(c) where the passing orbit in (r,9) coordinates are shown 
not to intersect the dashed line corresponding to the radius r tip , when ip = V^. The 
GC poloidal flux surface coordinate if) may be approximated by its orbit averaged value: 
(^orb = (' l P)orb(T :> (p, w , \ °~) ■ This idea has been given by P. Angelino et al. [4], who also 
suggest the estimate (ip) or b ~ 4>o — ^<j> ~ ( m sRo/ ( ls)o'V2w\ / l — XB H(1/_B — A), where 
the mass m s and the charge q s refer to the examined species, Rq and Bq are the major 
radius and the magnetic field magnitude at the magnetic axis, a = sgn(v\\/v) and the 
Heaviside function H(l/B — A) ensures that the square root is well defined. A more 
precise estimation will be given in (33) or (34) in Section (3.1). In [4] it is also suggested 
to slightly modify the equilibrium distribution function from a canonical Maxwellian to 
a biased canonical Maxwellian (where b means bulk and v t b = \jTb/mb): 

It is worth noting an inconvenience: having a density in tp means assigning a given 
number of GCs with a specific ip value that can be obtained from several V</>, w, A 
(and a). These QIs correspond to different orbits that span a wide portion of the GC 
configuration space. It becomes quite difficult to initialize that distribution function in 
a marker loading subroutine in a gyrokinetic code. Obviously this kind of problems are 
commonly addressed by ignoring the initial distribution function but possibly enhancing 
the number of markers to reduce the statistical fluctuation noise. 

A more serious problem consist in the reproduction of the experimental profiles as 
for the experimental density profiles n exp {jp) (for the temperature the analysis is even 
more complicated). Indeed, the fine property n exp = n e of the Maxwellian distribution 
(17) is definitively lost. This happens because ip, a spatial variable, has been substituted 
by ipQ that depends also on the GC velocity space variables (w and A). The same problem 
arises with the canonical Maxwellian distribution function because also mixes spatial 
coordinates with velocity coordinates. The properties of the Maxwellian distribution 
function are destroyed also considering (or ^o) as a spatial variable independent from 
the velocity space variables. Indeed the Jacobian necessary to achieve this independency 
will anyway break the characteristic Gaussian behavior of the integrand in the velocity 
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coordinates. Neverthless it seems very useful, as will be shown in the next section, to 
consider V<j, as an independent spatial coordinate! . 

Up to now the Maxwellian-like distribution function has been described. It is 
worthy to analyze other equilibrium (or steady state) distribution functions which are 
useful also for describing fast particles coming from fusion reactions or from external 
heating sources. The following model distributions are often used: the Slowing Down 
(SD) distribution function for fusion alpha particles 

where Ts is the Spitzer SD time [20] and w c = v 2 /2 is the critical energy [21, 22] and 
S a is the source term which corresponds to the density; the anisotropic SD distribution 
function for suprathermal ions coming from NNBI (Negative Neutral Beam Injection) 
heating [23] 

= T f j( A t\ H fr~:l e-^ 2 /^), (21) 
87r v / 2vrA(V') w 3 / 2 + vol' 2 

where £ = v \\/ v , w i is the beam energy, Sn{ip) is the source term and A (■?/>) gives the 
spread of the pitch angle distribution centered at £ ; the single pitch angle ICRH (ion 
cyclotron radiation heating) distribution function [24] for the minority population 

hcnni^X) - 27r2BovL(mm [—^j 5(\-X )e , (22) 



where 9b is the bounce angle near the tip of the banana orbit, Vt m = yT m /m m , T(z) is 
the Gamma function and n m (ip) is the density (m stands for minority); the modified 
biMaxwelian distribution function [25] 

3/2 ( 



f2Af(ip,w,\) = n h (ip) 



2ttT ± (V) 



exp < —rrihW 



\B res 1 1 — \B r 



(23) 



useful when the pressure tensor is diagonal but anisotropic {p» ^ p±) and where is 
the mass of the considered hot species, B res is a resonant magnetic field, nhijp) is the 
density and the temperatures T± and Ty can be deduced from the high energy limit 
respectively when \B res = 1 and when A = 0. 

The procedure analyzed before may be used also for the following distribution 
functions: the functional form is preserved and ip has to be substituted with its orbit 
averaged value (ip)orb (as done in [3] for the SD case). A similar prescription has 
to be used for the other evolving variables as £ substituted with (£) or b and so on. 
What is guaranteed with this ansatz is the dependency on QIs and the recovering of 
the commonly used distribution functions in the ZOW case, at least as regarding the 
dependence within vp. However this is not the only way to proceed. 

|| This is not a conceptual difficulty for those who works in gyrokinetics theory where also the GC 
spatial position is the difference of the particle spatial position minus a gyroradius that depends on the 
perpendicular velocity. 
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In this work it is proposed an alternative construction of the equilibrium distribution 
function which guarantees the above properties, but also the following ones: J- eq will 
behave according to a Boltzmann-like parametric distribution function and it will 
preserve useful integrability properties in the FOW case. How to build T eq will be shown 
in Section 4 whilst in the next section the orbits behavior will be described. Indeed, as 
mentioned before, J 7 eq (V ( p,w, A) describes a distribution of GC (projected) orbits. The 
analysis of orbit theory issues will serve to better understand the proposed equilibrium 
model distribution function. Moreover, in the next section the second condition (2), 
which is mostly ignored in the simpler theoretical models but useful in more realistic 
tokamak contexts, will be dealt with. 



3. Orbit theory fundamentals 



The unperturbed GC projected orbit is easily described with the following set of 
independent variables: tpjV^w and A (and eventually a). This system of reference is 
singular in the ZOW case, when = ip. However when the orbit behavior is considered 
the FOW effects must be taken into account. Using the definition of = ip + Fv\\/u c 
and the relation v 2 — 2w(l — A|S|), it is obtained: 



(n - v) 2 



2wF 2 (l - \\B\) 
IS- 



with F = m s F/q s . 



(24) 



Multiplyng both sides with B 2 , the unique positive solution of the second order equation 
in \B\ is denoted by B or i>: 

1/2 



B, 



w\F z 



orb 



1 + 



2(7^ -# 
w\ 2 F 2 



- 1 



(25) 



B orb is the intensity of the magnetic field magnitude \B\ seen from the GC along its 
orbital motion. It is worth noting that B orb depends only by tpjV^w and A which 
substantiates this choice of orbit coordinates. Once the (r, 9) map of the magnitude of 
the magnetic field \B\(r,9) is known, it is possible to describe the projected orbit in 
poloidal coordinates from the implicit relation: | B | (r, 9) = B^i^ir), V^, w,\). 

As an example, it is considered one of the most analyzed model [26] of tokamak 
plasma with nested circular flux surfaces with a Shafranov shift A(r) and a little 
inverse aspect ratio e, described with Shafranov coordinates (see Appendix A. 2), when 
A' = 0(e) (A.30, A.31): 

rF[l 



B = 



qR(R - A)(l 
F 
R 



0(e 2 )] F 
A'cosfl) 6 ' + R ( 



->• \B\ = 



1 + 



+ 0(s 3 ) 



(26) 



2q 2 (R - A) 2 

where e e = (Vip x V<f>)/\Vtp x V0| is the poloidal unit vector and q(r) is the safety 
factor. From (26) with _B orfe in place of \B\ and from R = R — A(r) + rcos9, it is 
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possible to express cos9 as a function of ^.Vs^w and A: 



cos 9 



rB, 



orb 



1 + 



+ 0(e s 



A-Rn 



(27) 



2g 2 (i? -A) 2 

Known the functions i[j(r), A(r), q(r), F(ip) and given P^, u> and A, (27) is the orbit 
expressed in poloidal coordinates (r, 9). 

Returning to the general |-B|(r, 9) case, it is possible to plot the orbit projection in 
the (r, 9) poloidal reference system if the QIs are assigned. It becomes easy to classify 
the GC orbits, also thanks to the chosen A coordinate. Indeed, the following relation is 
obtained equating B or b with \B\ and substituting w\ 2 with \ in (25): 

\ = XF 2 J f 2[V,-^(r)] 2 \ 1/2 

[V^-^r)Y\B\{r,9)\\ + X F* J 

No matter how complex it can be the map \B\(r,9), provided that the magnetic flux 
surfaces does exist, the A(r, 9) surface defined in (28) plays the same role of the potential 
energy in the classification of orbits in mechanics: the analysis which is based on the 
stationary points of the potential energy. 



A(r,0;7^,x). 



(28) 





(c) 



i.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0. 

rtip/a 



9 1 

r/a 



Figure 1. (a) Surface Polar Plot of A(r,0) defined in (28), at V$ and \ = w ^ 2 
fixed. The (projected) orbits are the level curves of A also recognizable in the Contour 
Plot just above the surface. As for example, the banana shape of the banana orbits 
is easily recognized. (b)|£| = versus r/a for the same case considered in the 

left figure. When |£| = the GC approximately reverse its motion, in which case the 
orbit is considered trapped, (c) The projected orbits plotted on the r/a — 6 plane are 
level curves of the A surface (on the left). It is possible to deduce the tip coordinates 
{Tup, Qtip) for each of the visualized trapped orbits from the intersections of those orbits 
with the dashed vertical line corresponding to £ = 0. 



Given a simple map \B\(r, 9) as (26) can be, and fixing and Xi ^ ne surface A is 
plotted in the polar (r, 9) reference system in Figure 1(a). Assigning a value to A means 
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cutting horizontally the A surface. The GC projected orbits are recognized as the level 
curves of A as clarified in Figure 1(c). In the depicted case it is assumed an up-down 
symmetry for simplicity. From the shown level curves of A it is possible to distinguish 
the lost orbits, if they intersect the radius r = r(ip a ) = a, where ip a is the poloidal 
flux of the magnetic field for the last nested magnetic flux surface (as for the separatrix 
magnetic flux surface when a divertor is present). In the depicted curves, it is possible 
to qualitatively distinguish the trapped orbit if they make a loop, from the others: the 
passing orbits, usually going from — n to ir, or viceversa. 

It can be also possible to quantify the above graphical method. For classifying 
orbits, the most common practice is to analyze the sign of the variable £ = v»/v = 
ay/1 — \Borb along the orbit: if £ changes its sign then it admits a zero and the GC 
orbit will be almost trapped, otherwise it will be almost passing. These differences are 
depicted on the right of Figure 1 comparing (b) with (c); Figure 1(b) shows the absolute 
value of £: 

The radius r tip and the angle 9 t i P are the coordinates corresponding to £ = at given 
V<p and X- When £ = then r t i P = r(ip) \^=-p tj> . Depending on A, it is possible to obtain 
or not, as the case may be, the intersection of the orbit with the r tip value (dashed line 
in Figure 1(c)). When the intersection happens, 9u p can be evaluated from the implicit 
relation: A = A(r tip , 9 tip ; V^, x)- While 9 tip depends on V^,x and A, it is important to 
emphasize that r tip depends only on V^, showing a degeneracy on A and w (or x)- In 
the Figure 1(c), it is also clear that the bounce coordinates (r b , 9b), defined as the values 
corresponding to 9' = 0, slightly differ to the tip coordinates (r tip ,9 tip )%. Although 
(fb,9b) ~ (r t ip,9 t i P ) is a good approximation, it is possible to be more precise. The 
classification of the orbit topology will be derived here directly from the geometry of 
the A surface. 

Generally, at most only five values of A can determine the whole classification of 
the orbits, as can be recognized in Figure 1(a): \ bm in and \bmax are respectively the 
minimum and the maximum value of A at the boundary, when r = a. \ sm in, Kmax are 
the local minimum and the local maximum of A, and they are commonly defined as the 
stagnation points, A c is the critical value corresponding to the saddle of A. Sometimes, 
the critical orbit (when A = A c ) is called the pinch-orbit. The coordinates (r c , 9 C ) are 
the solution of A c = A(r c , 9 C ; V<^, x) when A(r, 9) shows the saddle point. 

These particular values of A depend on V,p and on x- When an up-down symmetry 
is considered, they are displaced along the equator characterized by the abscissa 
x G (—a, a) . In Figure 2(a), the A section along the equator have been plotted to 
better visualize Xbmin, Kmax, X smin , \ sma x and A c for the same values of and x chosen 
for the plots in Figure 1. It can be seen how \ sm i n and A c lie on the High Field Side 
(HFS) when 9 = tt while \ sm ax lies on the Low Field Side (LFS) when = 0. 

% The difference between bounce and tip coordinates can be traced back to depend on the V|£?| drift. 
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(a)-l "°- 5 °- 5 x/a 1 (b)° 02 04 06 08 r/a 1 

Figure 2. (a) Equatorial section of A in Figure 1, in order to illustrate the A values 
used for the orbits classification: \bmim Xbmaxi X sm in, X sma xi A c , as described in the 
text.(b) (A, r) domain of variability subdivided in four zones: the loss orbits are 
enclosed in the white zone, the co(counter)-passing orbits for ions(electrons) are in 
the red(dark grey) zone, the trapped orbits are in the yellow(off-white) zone and the 
counter(co)-passing orbits for ions(electrons) are in the green(light grey) zone. 

The Figure 2(b) shows the same A section, once it has been folded around the 
x = axis. The area enclosed between the two branches A(r, 9 = 0), A(r, 9 = tt) and the 
vertical line r = a is the (A, r) domain of variability. The projected orbits are horizontal 
lines which connect the boundary of that domain. Now, it is possible to distinguish 
unambiguously the following four classes of orbits: (i) the loss orbits are those that touch 
the vertical line r = a, and are indicated in white. For the considered case, this happens 
for A < Xbmax, on the contrary, the confined orbits have A > Xb max (zones with colors). 
The confined orbits can be divided in (ii) trapped orbits if A > A c , indicated in yellow 
(off-white for b/w copy), or passing orbits if A < A c . The passing orbits are further 
subdivided between (iii) those with a radius r < r c , the counter-passing (co-passing) for 
ions (electrons), indicated in green (light grey for b/w copy), and (iv) those with a radius 
r > r c , the co-passing (counter-passing) for ions(electrons), indicated in red (dark grey 
for b/w copy). 

Figure 3(a) reproduce the same plot of Figure 2(b) together with the projected 
orbits plot of Figure 1(c), to better comprehend the relation between the classifications 
of orbits in the (A,r) domain and the shape of the corresponding orbits in the real 
space. The real space has been colored indicating the four zones: counter-passing(co- 
passing) for ions (electrons) in green (light grey), trapped orbits in yellow (off-white), co- 
passing(counter-passing) for ions(electrons) in red (dark grey) and loss orbits in white. 
Those zones are separated by the two branches of the critical orbit A = A c and by the 
confined boundary orbit A = Xbmax- In addition, to better appreciate the accuracy of 
the commonly used method based on the sign of £, it is also plotted the value of |£| 
(same plot of Figure 1(b)) and the dashed vertical line corresponds to r tip . 
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Figure 3. (A, r) domain of variability (as in Figure 2(b)) above the corresponding 
orbits plotted on the r/a — 8 plane as level curves of the A surface (as in Figure 1(c)), 
with the same color code used in Figure 2(b). In overlapping it is plotted |£| as in Figure 
1(b) and the r t i V value (dashed line). The employed V<p, X values are respectively: (a) 



V lj} = -0.5 Wb and X = 2 x 1CT 5 c 2 (b) 
and (c) = -0.5 Wb and X = 2 x 10~ 6 c 2 T" 



-1.3 Wb and X 



2x 10" 5 c 2 T" 



The plots that follow are realized in the same manner as for Figure 3(a), but A 
is now obtained with different values of V</> and x- m Figure 3(b), V</> is lowered. 
As a consequence the trapped orbits zone is shifted on the right (r c is increased as a 
consequence of taking ip' negative). In Figure 3(c), x have been lowered and V$ retained 
as in Figure 3(a), to show how the width of the orbits depends on x- the orbit width 
increases or decreases together with x- 

In Figure 4(a), x is kept as in Figure 3(b), but V</> is further lowered with the 
consequence that part of the trapped orbits are now lost: the confined zone becomes 
disconnected, separated by all the orbits with A c < A < Xbmax- 

When Vcj) is decreased below ip a , as depicted in Figure 4(b) and 4(c), there is place 
only for co-passing ions (or counter-passing electrons). The true reason is that A doesn't 
now show any stationary local maximum: X srna x doesn't exist in the plasma volume. In 
Figure 4(c), x is increased to show how r c can differ from r tip , indeed, in the present 
case, r tip doesn't exist because £ doesn't vanish. Whilst r Up depends only on V^, r c 
depends also on x- 

In Figure 5(a), x is kept as in Figure 3(a), but V$ is increased until \ S min is absent. 
The critical orbit has only one branch which is shown to separate the trapped orbits 
from the co-passing ions (counter-passing electrons). There is no place for counter- 
passing ions (co-passing electrons) because there aren't orbits with A < A c and r < r c . 
In the present case, the difference between r c and r tip is such that the common method 
of orbits classification, based on the sign of v\\ , cannot be applied without discrepancies. 

When V$ > ip{r) \ r =o= the A geometry changes further: \ smax is the only 
stationary point. In this case all the confined orbits are considered co(counter)-passing 
for ions (electrons). The Figures 5(b) and (c) differ only on the x value but they represent 
a similar case. 
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Figure 4. Same as Figure 3 to show the case where confined co-passing orbits for ions 
(counter-passing orbits for electrons) are not allowed. The employed V<j),x values are 
respectively: (a) = -1.73 Wb and x = 2 x 1CT 5 c 2 T~ 2 , (b) = -2.6 Wb and 
X = 2 x 1CT 5 c 2 T~ 2 and (c) = -2.6 Wb and X = 2 x 10~ 4 c 2 T~ 2 . 



(a) 



P«,= -0.14Wb (b) 
X = 2x10" 5 c 2 "T 2 



P* = 0.1Wb ( C ) 
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Figure 5. Same as Figure 3 and Figure 4 to show the case when confined counter- 
passing orbits for ions (co-passing orbits for electrons) are not allowed. The employed 
"P^xvalues are respectively: (a) V4, = -0.14 Wb and x — 2 x 1CT 5 c 2 T~ 2 , (b) 
V4, = 0.1 Wb and X = 2 x lO" 5 c 2 T~ 2 and (c) = 0.1 Wb and X = 2 X lO" 3 c 2 T~ 2 . 



3.1. Orbit average, characteristic frequencies and second invariant 

In the orbit coordinates i(},V^,w and A (and eventually a) it is possible to express the 
orbit average. The orbit average of a quantity A is defined as a time average along the 
path C which is the closed projected orbit (only the projection of the orbit is always 
closed): 



J c Adt _ f c Ad6/6 _ J c Adr/f _ J c A u r 
{ )oTb Icdt Ic d9/9 fcdr/r ^ 



It is common practice to substitute one of the relations (3) or (4) to evaluate the integral. 
An alternative procedure is here adopted which uses the relations (16) expressing the 
QIs. Thus, the QI average of A is defined as follows: 

(A) = f A QrHjVt - V^)5(w - w)6(X - A) d 3 id 3 v 

{ )QI j QrmPt - V^)6(w - w)5(X - A) d 3 xd 3 v ' 1 j 
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where the integral is evaluated on the whole phase space when the values of V^, w and 
A uniquely determine the orbit. When there is a degeneracy, as in the case of the co- 
and the counter-passing orbits, the phase space must be divided in disjoint subspaces 
where the uniqueness of the orbit is recovered. These subspaces are distinguished by an 
identifier a. 

In (31) the volume element comes to be 

- wlB }^ d^dedwdXd^ = ^^ygd^dwdx^dnf 

V\v\\\ W N 00^1 

with y/g and dgV^ computed at 9 — 9(il),V<t>,w, A), when \B\ = B or t,; 7 is the gyrophase 
which is an ignorable coordinate as well as the toroidal angle 0. In order to demonstrate 
the equivalence of (30) with (31), the first equivalence of (11) is used in (32). Finally, 
the obtained volume element is used into (31), deducing the equivalence of the QI 
average with the orbit average (30). As an application of (31), here it is the value {ip) or b 
requested for the biased canonical Maxwellian distribution function (19) described in 
Section (2.2): 

where ip m in and tp max are respectively the minimum and the maximum poloidal magnetic 
flux value reached by the orbit. Replacing ip = V<f, — Fv\\/u> c : 

ii\ =v _ rn,J$Z:FJg/WQd B V+)dil> 

K l0Tb * VsSt-B^g/WQvwdeV^d^ 1 1 

where the approximation {ip) or b ~ ipo used in [4], can be deduced when \B\ = B orb ~ B . 

Another application of the above expression for the time integration is on the 
computation of the characteristic frequencies, the inverse of the bounce time r b of the 
projected orbits: 

n = f dt= [ 5{% - V+Mw - w)8(X - A) (35) 

JC J 47T SzlW 

finding out: 

n = _ 2 f^ max B orbV9 d ^ = 2 jqs[_ f r ™* VgBjbdr 

U™ n ^'Q\v\\d e V^\ wm s Jr min QF(2 - XBor^ldelBW' 

where the relation v^doP^ = —wFB~ 2 [2 — X\B\]d e \B\ has been used. In the above 
integral y/g, Q and dg\B\ have to be computed for 6 = 0(il),V<i>,w, A), when \B\ = B or i>. 
Moreover, r min {%l) min ) and r max (ip max ) are respectively the minimum and the maximum 
radius (poloidal magnetic flux) value reached by the orbit. 

It is worth noting that the bounce frequency and the transit frequency are the same 
frequency 2n/T b (V ( i ) , w , A) computed respectively for A > A c and for A < A c + . One 
can show how the critical frequency 2ir /^{V^, w, A c ) goes to zero, if A c is allowed once 

+ The difference between the co-passing and the counter-passing orbits is on the integration range of 
the radial coordinate. 
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V(j, and w are given. The limiting case of the orbit behavior in close analogy to the 
pendulum, analyzed in details by Brizard A. J. et al.[27], can be recovered, as shown by 
Chiu S. C. et al.[15]. 

A quick remark is concerning the second invariant J, often used as COM instead of 
w, with the property: d w J = r b . It is clear from (35) how it can be expressed in terms 
of V</,, w and A: 



3.2. How to select only confined Guiding Centers 

The reasoning for considering only confined orbits is usually fairly structured [11, 12]. 
A simple condition that may be implemented in the codes, as done in [28], is analyzed. 
For simplicity, it is assumed that V\B\ ^ 0: for each r (each tp) there exists only 
one maximum HFS(r) and one minimum LFS(r) of \B\(r,6) and these are the only 
stationary points when de\B\ =0 (the magnetic flux surface is convex). In this way it is 
possible to define an equator as the locus where \B\ = HFS(r) on the High Field Side, 
passing through \B\ = Bq at r = 0, and where \B\ = LFS{r) on the Low Field Side. 
The value of 6 can be fixed to 9 = in the LFS direction and to 6 = n in the opposite 
(HFS) direction. These rules become evident for the up-down symmetric plasma where 
the plane of symmetry is the equatorial plane. 

As a starting point it is analyzed the range of variability of the following quantities: 
A G [0, 1/B min \ : \B\ G [B min: B max ] where B 

max — disx^j) \B\ and B m i n — min( rj g) \B\. 
Moreover, ip G [V'ojO] corresponds to r G [0,a], being ip(r) \ r =o= and ip(r) \ r=a — ip a < 
0. Using the expression (25) of \B\ in orbit coordinates the condition B min < B^b ^ 
B max must be always verified. If it happens that B m i n < B or b(ip a , V^>, A, w) < B max then 
the orbit can intersect the boundary of the plasma volume, the r = a surface, and the 
GC will be lost. The condition for the confined GC will be the complementary one, 
when one of the two situations occurs: 



The critical case when the orbit is tangent to the surface r = a has been considered as 
if the GC is lost. It is possible to rewrite the above condition as follows: 



because the maximum and the minimum value of A at the boundary surface r = a are 
respectively: 




(37) 



B = {(B orb (ip a , V<f„ A, w) < B min ) or {B orb (ip a , A, w) > B max )}. (38) 



B = {(A > XbmaxCPtf^x)) or (A < Afe min (P ,x))}, 



(39) 




(40) 



and 




(41) 
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where F a = F(ip a ). The condition (38) is easily visualized as satisfied in Figures 3(a) 
and (c), in Figure 4(b) and for the Figures 5(a),(b) and (c). 

In the favorable case in which a particular set of (V,p,w,X) selects only a single 
orbit, the condition (38) or (39) provides that the orbit will be confined. The problem 
arises when the same set of orbit coordinates represents two orbits. If these orbits are 
both confined or both loss (hardly), the above condition will be respectively satisfied, 
as for Figures 3(a) and (c), or broken. However, when one of the orbits is confined 
and the other is not, the condition (38) (or (39)) is not anymore sufficient; as it is in 
the cases of Figure 3(b), Figures 4(a) and (c). Fortunately, it is easy to discriminate 
which of the two orbits is confined. Indeed, it is possible to have two orbits when the 
surface A admits one saddle point at A c = A(r c , 9 C ) and A < A c , when the orbits are 
passing. While the critical A c divides the orbit in passing or trapped, the critical radius 
r c divides the passing orbits in two disjoint families of orbits: for r < r c the orbits are 
counter-passing for ions and co-passing for electrons, viceversa for r > r c . 

It is obvious that the orbit will be confined if the condition 

A = {3(r c , A c ) : (r < r c (^, *)) and (A < A c (^, *))}, (42) 

will occur, regardless (38) or (39). 

In summary, the condition for taking into account only confined GC is: 

| 1, if AUB 

0~ confined 1 . (43) 

I 0, otherwise. 

The criterion (43) is easy to implement in a gyrokinetic code thanks to a Metropolis 
algorithm which is explained below. Gyrokinetic codes provide an initial value to the 
GC coordinates taking care of the uncorrelated initialization between a couple of GCs. 
For example it is used a particular set of coordinates: Z = (r,0,(f>,v\\, [/). Then it is 
necessary to define the subset of confined GC depending on V</>.w, A and a: 

Qconfined {'Pcfi-U], A, 0~ | 5 C onfined I}- (44) 

Initialization starts with the assignment of an array of values Z — Z\. From these values 
it is possible to calculate V<f>i = V^[Z\),w\ = w(Zi),\i = \i(Z~i) and o\ = o(Z\). If 
Qi — {^ > <t>ii w ii Ai, (7i) G Qconfined then the value is taken and the initialization proceeds 
with the next Z2. Otherwise the set of values is rejected and the code returned to 
reinitialize Z 1 with another set of values. Up to load all the GCs. 

4. A new parametric distribution function for gyrokinetic equilibria 

Three guidelines are required to proceed on building the equilibrium distribution 
function T eq : (a) it behaves as a Boltzmann distribution function, T eq oc exp — £/T; 
(b) £ is function of QIs; (c) T eq is mathematically tractable. The first request, although 
only a guideline, concerns a property that suggests to consider T eq something more 
than a simple fitting model distribution function. Indeed, the Boltzmann behavior 
commonly expresses a physical property: it is expected when the two-point correlations 
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between particles can be neglected [29]. The second issue is mandatory once assumed 
a Boltzmann behavior, because (1): T eq has to depend on QIs to be constant in time. 
The third request is speculative, but very useful. 

In what follows, it is shown how it is possible to arrive at an expression of the 
energy as a function of COMs. If the distribution of particles is described by a local 
Maxwellian, many of these will have |£| ~ 1/3 for the reached isotropy. The same for 
the SD distribution function. When it is considered a minority species from ICRH there 
will be an anisotropy such that it is easier to find a particle with |£| < 1/3, or else, when 
they are considered energetic particles coming from a beam there will be an opposite 
anisotropy such that |£| ~ 1 is favored. The ranges |£|< 1/3 and |£| ~ 1 are considered 
with more care. 

Here it is convenient to switch on using particle coordinates, although the symbols 
used for GC coordinates are preserved. As a starting point it is considered the kinetic 
energy of the single particle: 

with sin a = C, = v\\/v and |e| < 1 ( v\\ ~ when \B\ ~ B<p as in tokamaks so that 
e = (arcsin^/v)/(arcsin v\\/v) — 1 is little enough). When £ is properly small to allow 
the expansion a ~ £ + £ 3 /6 then the following term can be expanded: 

cos 2 [a(l + e)] ~ 1 - a 2 (l + e) 2 + \a\l + e) 4 



£ 2 (l + e) 2 + C(2 + e)(l + e) 2 = 



= 1 3(1 + 6) 2 | e(2 + e )(l + e ) 2 



3 

i _l Aft _l A2 r 

2 



2e(2 + e) 



4e(2 + e) 3 

= l-a + b (£ 2 + c -l) 2 , 

where a = 3(1 + e) 2 /[4e(2 + e)],b = e(2 + e)(l + e) 2 /3 and c = 1 - 3/[2e(2 + e)], for 
convenience. Then, substituting 

vj {Ru^_ (L z /g s - ^ 

2 " 2B? 2{m s R/q s y ' 1 ° J 

(45) can be rewritten as follows: 

How it is evident, this relation follows from the balance between the first and the third 
term, the sum of which gives an almost zero contribution. Inserting a new parameter 
k > max(l,a ) then adding and subtracting kv 2 /2, it is obtained 

T~TO? + £<i-«> + £[«-«. + MA|*l -*)*]. (48) 
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Finally, putting on the LHS the second term it is possible to rescale the energy, now 
unbalanced between the component and the one orthogonal to it: 



^ (L z /q s - jj) 2 w(k - op) 
2 2n(m K R/a f! ) 2 k 



1 + 



b \B\(\-c /\B\f 



(49) 



2K{m s R/q s ) 2 k k — ao 

A clear functional dependency on (w, A, Lz) of the energy £ ~ m s v 2 /2 has been found 
when a ~ £ + £ 3 /6 is a good approximation. The case in which ( ~ 1 or a ~ tt/2 is 
studied on the following. A slightly different expression respect to (45) is assumed: 

,2 



V 

~2 



v l + v l 
2 2 



r 7T 



7I\ 



sin 



[o+(«- o)(l-t)] 



(50) 



■2 v 2' 

with the same condition |e| < 1. Expanding the sine up to the forth order on (a — vr/2) 
and preserving terms up to the (1 — £ 2 ) 2 order, the same steps as before are followed, 
to arrive at 



v 
~2 



(L z /q s -ijf 



W{K - aij - ci/ |^|; 

2(l + /t)(m s i?/g s ) 2 (1 + k) [ « - a x J' 1 J 

with ai = 3(1 - e) 2 /[4e(2 - e)]A = e(2 - e)(l - e) 2 /3, ci = -3/[2e(2 - e)] and 
/t > max(l,ai). A similar result is also obtained for £ ~ — 1 or a ~ —n/2. The 
above relations are all expressed in particle coordinates. (49) and (51) will change 
once expressed in GC coordinates. Without entering into the details of the GC 
transformation, it can be assumed that the given expressions are preserved up to the 
lowest order if transformed in GC coordinates. Five parameters (V^o, Ao, Ap , T w , A\) 
that depend on the GC position, are introduced. Lz/q s is substituted with and £ 
is divided by T a reference constant temperature per unitary mass,. Finally, both (49) 
and (51) become (after the GC transformation, in GC coordinates): 



6i|S|(A-ci/|S| 



£ 
T 



in - V, 



<t>o) 



A Ps 



w 

+ 

T 



1 ± 



(A - Ao)- 
A! 



(52) 



with the squared parenthesis defined positively, such as A% , , T w and A 2 . 



A GC Boltzmann-like equilibrium distribution function can now be written down 
using (52) in the exponent. To fulfill the first equilibrium condition (1), the functional 
dependence on the QIs is maintained, whilst the not constant quantities are replaced 
with constant parameters to be determined afterwards. In this way a parametric 
equilibrium distribution function is obtained. It leaves to the parameters the task to 
capture the collective character of how the GCs are distribuited. Indeed, regardless 
of the plasma operational configuration and of the particle species (fusion products, 
thermal bulk, energetic particles from ICRH and NNBI), the single particle energy is 
always given in the form £ = m s v 2 /2 + q s A . This consideration leads to a unique 
Boltzmann-like distribution function to be valid in any context. However, ensemble 
phenomena arise: a collection of these GCs is distributed in a non-uniform fashion; 
e.g., concentrated in the hot core plasma as for fusion alpha, concentrated off-axis as for 
ICRH, or having a narrow pitch angle distribution with dominant v» as for NNBI. These 
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Figure 6. (a) Flux surface averaged density profile n in a.u. versus r/a for 
the canonical Maxwellian case computed from the dist. func. (54) with T w = 
5.0 x 1(T 6 c 2 , 7> = 0.0 Wb, A 2 P = 0.1 Wb 2 . (b) Surface polar plot plus contour 
plot of the density n in a.u. versus (r, 9) computed from the same dist. func. in (a). 
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Figure 7. (a) Flux surface averaged density n profile in a.u. versus r/a for the ICRH 
case computed from the dist. func. (53) with a = 1.25, T w = 0.00025 c 2 , Ao = 0.13 
T"\ A\ = 0.00001 T- 2 , -P 0O = -1-0 Wb, = 0.05 Wb 2 . (b) Surface polar plot 
plus contour plot of the density n in a. u. versus (r, 9) computed from the same dist. 
func. in (a). A lot of GCs are positioned near the tips of the banana orbits, forming 
two horns. 



peculiarities must be taken into account, firstly, separating the perpendicular from the 
parallel anisotropy component in the kinetic energy (as done above introducing the 
unbalancing k parameter for the particle kinetic energy expression), then by choosing 
the appropriate parameters. 

As an example, the following distribution function is taken into consideration: 



eq 



/27TW 3 / 2 

where J\f,T w ,a,V< 



w 



exp 



V A -V, 



©0 



exp 



w 

in 



1 + 



(53) 



Ap^, Aq and Aa are constant parameters. 



The canonical Maxwellian with constant Temperature (and a Gaussian behavior in 
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V4) is obtained for A\ — > oo and for a = 3/2, as shown in Figure 6: 



lim F eq (a 



3/2) 



Af 



2^ /2< 



cxp 



cxp - 



w 

-L it; 



(54) 



-w/T w 



A local Maxwellian in the ZOW limit is obtained: 



lim F zow {a = 3/2) 



2ttT, 



3/2 



cxp 



'ip-V, 



<t>o 



cxp - 



w 

J- in 



(55) 



With the case A A ->■ 0, Af = A"A A y^0,a = 5/4, (53) becomes 



lim F eo (a 



5/4) 



AT 



V2T, 



3/2 



T \ 3 /4 



cxp 



(56) 



exp - 



w 



T 



(J(A-Ao), 



which is in accordance with (22) in the ZOW limit: 



lim -F zcw (« = 5/4) 



AT 



3/2 



T 



w 



3/4 



exp 



(57) 



exp - 



w 



-*- ii 



S(X-Xo). 



Especially, when the FOW effects are considered, it is observed the expected minority 
phenomenology in the presence of an ICRH antenna. For example in Figures 7 and 
8, a higher concentration of minority near the banana orbit tips positioned along the 
resonant value of the magnetic field magnitude B res = Aq 1 is shown [30]. 

Finally, the modified biMaxwellian distribution function (23) can be partially 
reproduced. Even when FOW effetcs are taken into account such similarities are 
observed when (53) is rewritten as 



3/2 



-m h w 



F eo (ct = 3/2) = n h (Vd),X) exp 

qK 11 V ; (2vrT ± (A)) 3 / 2 l 

where the following definition has been adopted: 

m h A\T w 



\]3 res 1 \B res 

tax) W) 



(58) 



and 



f±(X) 



A\ + A 2 + Ag 

m h A\T w 
AJ + A'-Ag' 

2%AfA* 



rhf) 
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The similarity with (23) is only when A < Ao = B rf } s . On the contrary, the following 
form of the equilibrium distribution function should be considered to obtain (23) when 



A > A = B: 



l. 



eq 



exp 



V, 



00 





w 


exp < 





A-Ar 



(59) 



with the minus sign in the squared parenthesis of the last exponent. This case can be 
obtainable from (52) too. Even though interesting, in the following the analysis will be 
focalized on the equilibrium distribution function (53). 



4-1. Useful mathematical properties of the proposed distribution function 

The equilibrium distribution function (53) is a totally parametric function particularly 
indicated for the differential calculus, as can be the application of a differential 
collisional operator. It worths the trouble to emphasize the integral properties of 
the obtained distribution function mainly used for the velocity momenta computation. 
Furthermore it is important to justify the choice of the normalization factor proportional 
to (w /T w ) a / \ / 2ttw 3 ) in front of the Boltzmann-like exponential law in (53). Going back 
to the well known variables (ip, 9, w, v \\): 
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(53) can be written as follows: 
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exp 



exp 
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In the above expression the —v^/w term is crucial (the sign justifies the choice of the 
plus instead of the minus sign in (53) respect to (59)). When the only v\\ dependence 



in F eq is considered: 



Feq{v\\ 



M exp -Jv\\ 



2 2 

m vf\ 



9^1 
4! 



(62) 
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r/a n(a.u.) r/a n(a.u.) 




Figure 8. Contour plots of the spatial density in r/a — 8 plane for the ICRH case 
computed from the dist. func. (53) with: (a) a = 1.25, T w = 0.0002 c 2 , A = 0.13 
T"\ A 2 = 0.00001 T~ 2 , 7> = -0.85 Wb, A% = 0.02 Wb 2 ; (b) a = 1.25, T w = 
0.0002 c 2 , A = 0.13 T-\ A 2 = 0.000008 T~ 2 , V^o = -1-0 Wb, A 2 = 0.01 Wb 2 . 
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The exponent is very similar to the action in the usually called "A0 4 " field theory [31], 
with an interaction term " J ■ 0" where J is the interaction current (when m 2 < and 
J = the classic double- well potential is recognized). Thanks to this coincidence, 
it should be possible to borrow for our use some techniques used in quantum (or 
condensed matter) field theory. The following formal identity can be useful when 
r = (AwT w B 2 A 2 ) -1 is sufficiently small and q ^ (or m 2 ^ 0): 
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The chosen factor in front of the exponent in (53) is justified when one explicits 
the only w dependence. F eq is rewritten as follows: 



F eq (w) 



Af ( w 
2W/2 KJ], 



/ aw b 

exp 



with 



and 



Af = Af exp 



exp -Jv\\ - 



1 2 > 
m v» 



(66) 



(67) 



(1-A |S|) 2 + S 2 A 



T W B 2 A 2 



(68) 



2T W B 2 A 2 X 



The function in (66) is recognized to be the statistical Inverse Gaussian (IG) distribution 
when a — 0, otherwise it is proportional to the Generalized Inverse Gaussian (GIG) 
distribution. The properties of IG and GIG pdfs (probability density functions) are well 
known [32, 33, 34]. For example, the definition of GIG pdf itself is used to show the 
identity: 



00 w n dw 
\/2nw 3 / 2 
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-*- tn 



exp 



aw 
'^2 



2K p (Va6) 



27iT«(a/b)p/ 2 



p = n + a -1/2, (69) 



where K p is the modified Bessel function of the second kind. 

These and other mathematical properties will not be used here. It was important 
to emphasize the powerful possibilities offered by (53) when a velocity integration has 
to be performed. However, it should be said that the determination of the velocity mo- 
ments can be obtained in an elementary way only when A A goes to infinity, otherwise 
the use of incomplete Bessel functions is required. 



4-2. Energy boundary behavior and confined Guiding Centers condition 

The distribution function (53) gives rise to a logarithmic divergency in the low energy 
limit (infrared divergency) for a = 0. A simple regularization scheme replacing w 3 ^ 2 
with w 3 / 2 + w^ 2 is applied to retrieve this interesting case. Looking at (9) and at (21), 
the parameter wq will correspond to the critical energy w c [21, 22]. Conversely, the 
ultraviolet limit (w — > 00) does not commonly apply because there is an upper limit on 
the usable amount of energy. In this way, it is appropriate multiplying the just obtained 
distribution function for a step function H(w\ —w) (or for some more smoothed sigmoid 
function). The cut on energy is given at W\. 

At the same time the condition selecting only confined orbits is symbolically 
represented with the factor 5 con fi ne d, defined in (43). 
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(a) 




Figure 9. (a) Surface polar plot plus contour plot of the density n versus (r, 9) in a.u. 
for the SD case computed from the dist. func. (71) with V^q — 0.0 Wb, A 2 ^ = 0.05 
Wb 2 . (b) Velocity density distribution d v in a.u. versus (w,v^) computed from the 
same dist. func. of (a). The ridge of the velocity density profile is the standard u>~ 3 / 2 
behavior. 



Finally the following expression for the equilibrium distribution function is 
obtained: 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.03 -0.02 -0.01 0.01 0.02 0.03 

(a) r/a (b) vn(n.u.) 



Figure 10. (a) Flux surface averaged density n profile in a.u. versus r/a fot the NNBI 
case computed from the dist. func. (73) with A = 0.08 T _1 , T W A\ = 5.0 x 10~ 8 c 2 
T" 2 , V^o = 3.0 Wb, A 2 ^ = 0.15 Wb 2 . (b) Contour plot of the velocity density 
distribution d v in a.u. versus (w,v») computed from the dist. func. of (a). 

It is possible to find some set of parameters useful to describe a slowing down of 
the energy. For T w — > oo and for A x ^ 0,a = 0, the SD-like distribution function (20) 
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is obtained with a GC anisotropy caused by having a (in place of ip) dependency: 



lim J- eq (a = 0) 



T w ^oo 



27r[w 3 /2 + w 3/2- 



cxp 



V, 
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H(iui - w)5 confined . (71) 



Figure 9(a) shows the spatial density in polar coordinates computed from (71). Figure 
9(b) shows the velocity density distribution in the (w,v\\) coordinates for the same 
distribution function. (71) becomes a SD in the ZOW limit: 
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t„^oo — V2n[w 3 / 2 + w 3 /2 ' 

Even more possibilities arise when T w — > oo and A A — >■ such as T W A 2 X ^ is finite. 
When a = 0, it is possible to obtain a distribution function that allows the modeling of 
ions heated by a NNBI: 

exp 
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(73) 



The "perpendicular anisotropy" term exp[— iy(A — Ao) 2 /(v / 7^Aa) 2 ], together with the 
"parallel anisotropy" term exp[— (V^ — V<f>o) 2 /A 2 P J, determines the desired anisotropy 
in the presence of a NNBI source. The A-term determines how many passing orbits 
are loaded in respect to the unwanted trapped orbits: Ao must be less than A c . The 
P^-term tunes the balance between co- and counter- passing orbits. Figure 10(b) shows 
the velocity distribution as function of (w,v\\) in arbitrary units: the energy follows the 
behavior w ~ v f\/2 which indicates that most of the GCs are passing (in the same plot 
it is shown an imbalance for the co-passing orbits). Figure 10(a) shows the flux surface 
average of the density even obtained from (73). It is worth noting that the anisotropy 
is mainly due to A in place of £ with respect to (21). A similar choice has already been 
proposed by [35] where the Figure 1 therein confirms the result in this article of having 
a pitch angle width A A that depends on energy: 



exp 
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5. Qualitative description of the equilibrium distribution 

The factors present in (70) are here described in a qualitative way. T eq is rewritten here 
for clearness: 
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The following terms can be distinguished. 

The Maxwellian term exp(— w/T w ) is well known, when Ax — > oo, Ap^ — > oo and 
w\ — > oo then T w becomes the temperature (per unit mass) expressing the decay rate 
of the energy. For wq = and a = 3/2 the Maxwell distribution function is retrived. 

The term (2tt)- 1 / 2 A![w 3/2 + w^ 2 ]' 1 indicates the SD of the energy. This behavior 
is widely known and occurs, for an appropriate range of energies, when modeling the 
SD distribution function for fusion alpha particles and also for the NBI distribution 
function for supra-thermal ions (in which cases w^ 2 ~ (ni/n e )v^/2 3 ^ 2 with v c the critical 
velocity [21, 22]). When the important range of energies is much higher then w , this 
can be ignored leaving the standard factor Af / V2ttuP of the IG pdf. M is the overall 
normalization constant which cannot be given explicitly because it is hard to estimate 
the integration of T eq on the entire configuration space despite the useful mathematical 
properties previously shown. 

The term (w/T w ) a is a power law used to take into account mainly the low energy 
behavior. It is fundamental to evaluate the content of energy that decreases when a 
increases. When a = and T w — > oo,Ap — > oo the SD distribution function (20) is 
retrieved. When a = 5/4 and A A — > 0, similarities with (22) are found. When a is 
an integer and wo = then it is possible to express the integrals on w with analycal 
combination of error functions. 

The Heaviside step function H(u>i — w) is used when the particles described are 
created with a given energy: w±, e.g. wi ~ 3.52 MeV for the alpha particles. In the 
presence of a beam, this term derives from the monochromatic source approximation 
(oc 5(w — wi)). When the beam cannot be considered monochromatic, as for source 
oc exp[— (w — u>i) 2 /A^], then it would be better using erfc[(u> — Wx)/A w ] instead of 
H{w\ — w). 

The term exp[— w(\ — A ) 2 / (T w A 2 x )j is a completely new term. The presence 
of the energy, together with the generalized pitch angle, prevents the factorization 
of the equilibrium distribution function as oc f(V < p)g(w)h(X). It is worth noticing 
the fundamental exception of the single pitch angle case, retrieved as the limit of 
\Jw/ (T w A"yj exp[— w(\ — Ao) 2 /(T w A 2 % )] for ^JTwA\/w — > 0. The relevance of A on 
classifying the topologies of the orbits has been pointed out above. This term can also 
counts the number of trapped orbits respect to the passing one. If there is an interest 
in studying the behavior of the passing orbits only, then A ~ has to be properly set 
(deeply passing orbits limit is when A — > 0). Moreover, to study only trapped orbits, 
one has to properly set Ao ~ \/B res for a given resonant intensity magnetic field. In this 
case a large amount of GCs is deposited into the region with B ~ B res . For example, 
this can be the case for simulating the minority (m m , q m ) in the ICRH scheme, for which 

UlCRH = Q.mBres/ m m- 

The exp[— (V^ — V^o) 2 / A 2 P } term is noteworthy and clear. V<p represents the 
projection of the orbit at given w and A (and a). This means orbits are distributed 
in the most simple way, that is with a Gaussian around a mean value V^o- In the SOW 
case V<f, ~ ip oc r 2 in the proximity of the magnetic axis. The tail of this distribution 
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will go rapidly to zero as exp — (r/A r ) 4 when V^o ~ 0. On the contrary, when V^o is 
taken outside the range of the allowed values of ip, the presence of only passing orbits 
is facilitated (see Figures 4 (b) and (c), or 5 (b) and (c)). In this case, it will result an 
imbalance on v\\ due to J in (63), such as in the case of studied population coming from 
a (N)NBI heating. 

The last S con fi ne d term should be used mainly when there are many orbits with a 
large width that can be lost. This term can be implemented numerically thanks to (43) 
as described in detail in Section 3. 

6. Conclusions and perspetives 

This work addressed the problem of the equilibrium distribution function in the 
gyrokinetic theory. It has been defined an equilibrium distribution function of GCs 
fulfilling the following conditions: (1) it must depend only on invariants of motion and 
(2) GCs must remain confined for suitably long time. The chosen set of invariants 
is (V<f,iWi\). These invariants have been called Quasi Invariants (QIs) as clarified 
in Section 2. The Section (2.1) emphasizes the connection of the expression of 
= ip + Fv\\/oj c with the expression of the drift velocity v D (5), thanks to the 
symmetry. Section 2.2 summarizes the way how address the currently studied 
problem concerning the equilibrium in gyrokinetic simulations in view of emphasizing 
the alternative approach used in this article. Some results of the orbit theory have 
been recoverd in Section 3 to introduce the reader to the (ip, V<f>, w, A) orbit coordinates 
used for toroidal symmetric plasma. Some results may be considered more accurate, 
for example a clear visualization of the orbits and their classification due to the surface 
A defined in (28). In addition it has been proposed a new way to compute the orbit 
average (31) that allows to write the formal expression of {ip) rb in (33), as well as the 
characteristic orbit frequency from the bounce time (36). The conditions to discern 
whether the orbits are confined or loss have been determined ((39), (42) and (43)). In 
Section 4, the equilibrium distribution function has been constructed in parametric 
form with the following three guidelines: (a) a Boltzmann-like distribution function, 
T eq oc exp— £/T; (b) £ has to be expressed as function of QIs; (c) T eq has to be 
mathematically tractable such as being used in integro-differential calculus. T eq has 
shown some asymptotic behaviors that are typical of some of the most used distribution 
functions in the gyrokinetic theory for tokamak plasma namely the local and canonical 
Maxwellian and the Slowing Down distribution function, respectively (55), (54) and 
(72). Moreover, the obtained J- eq shows analytic similarities with other distribution 
functions such as the single pitch angle, the anisotropic Slowing Down and the modified 
biMaxwellian distribution functions, respectively in (71), (73) and (58). It is given 
an explanation for the good comparison of the behaviors wanted from external sources, 
thanks to the possibility of selecting the kind of orbits which are mostly loaded: trapped 
orbit for the ICRH minority distribution function or passing orbits for the suprathermal 
ions from NNBI source. This can be deduced from the Figure 10(b) showing the 
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velocity density distribution for the anistopic SD distribution function and indicating the 
suppression of trapped particles respect to the passing one. On the contrary, the banana 
shape for the minority density in the ICRH case in Figures 7(b) and 8, is clearly due to 
the prevalence of trapped orbits. The Figures 8 show the density contour plots computed 
for the ICRH case, and they are surprisingly (because it is machine independent) very 
similar to the ones reported in [30]. In Section 4.1, the mathematical properties of T eq 
are partially described. These properties arise from the functional behavior on w and 
on v || that corresponds to the well known cases encountered in statistic analysis and 
in quantum (or condensed matter) field theory. A summary of the various factors that 
constitute the T eq has been qualitatively provided in Section 5. 

T eq can now be used to fit experimental profiles and it could provide a useful tool for 
experimental and numerical data analysis. The proposed model distribution function 
can be easily implemented in gyrokinetic codes because it is basically an analytical 
function. In Section (3.2) a method is proposed to implement the condition to avoid loss 
orbit in the loading subroutine of a gyrokinetic code. This model distribution function 
has already been used to simulate plasma in the presence of external heating sources, as 
demonstrated in [36] for the ICRH case relating to FAST [37] plasma conditions, using 
the HMGC code [38]. 

The importance of having a functional fully parametric form has not been shown 
here. However, the analytical advantages of applying differential operators to it (e.g. the 
collision operator) are evident as well as the possibility to perform parametric studies. 
For example, it would be possible to relate all the results occurring from a gyrokinetic 
simulation to the values of the six parameters: a, T w , V^o, Ap^, \ and A\. This is what 
it is commonly required to carry on the benchmarking of codes. 

A final point concerns the modeling of the distribution functions out of the 
equilibrium. Considering the simple case of the evolution of a plasma population 
distribution function as transiting between close equilibrium states, this can be easily 
addressed giving a time dependency to the parameters of T eq . While it may seem 
premature to consider a possible use of T eq in transport models, it seems equally clear 
that it can be applied to the integrated simulation of plasma scenarios [39]. The dialogue 
between the various codes would be highly optimized if the passage of information occurs 
through the relevant T eq parameters. 
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Appendix A. Magnetic flux and Shafranov coordinate system 

In plasma theory and modeling a very useful and often adopted representation of the 
magnetic field is the magnetic flux representation. 

The Shafranov coordinates are particularly used in the context of plasmas with 
circular poloidal flux section geometry. Once the Grad- Shafranov equilibrium equation 
[41, 42] is imposed to the system, this coordinate system leads to the s-a model [43]. 

In this appendix the magnetic field is briefly described in terms of flux coordinates 
before moving to Shafranov coordinates. Moreover, it is shown the correspondence 
between the two aforementioned representations. 

Appendix A.l. Magnetic field flux coordinates representation 
In flux coordinates B is simply written as follows: 

B = Vip x V0 + ^-V<pt x W, (A.l) 

where is the toroidal coordinate, i? is the poloidal angle in flux coordinate, —2tii/j = <p p 
is the poloidal magnetic flux and <p t is the toroidal magnetic flux. Thanks to V ■ B = 
and the Gauss theorem, the fluxes can be expressed as: 

1 " B-V§d 3 x (A.2) 



p 2n Jn(ti>) 

<t>t = 7r[ B-V<j)d 3 x, (A.3) 

where is the plasma volume enclosed into the magnetic flux surface ip. The toroidal 
coordinate system (r, 0), where the radial coordinate r labels the magnetic flux surface 
r = r(ip), is known as the flux coordinate system if 

d$ = d(j)/q(r) (A.4) 

along the magnetic field line. (A.4) is known as the straight-line property of the magnetic 
field line. In flux coordinates the volume element d 3 x becomes d 3 x = y/gfi ux drdi9d(j) 
where the Jacobian is 

Jg]^ x = (Vr x W • V0)- 1 = -tf/B ■ W, (A.5) 

where the prime indicates the radial derivative. 

In (A.4) q(r) is the safety factor which can also be expressed as: 
( s B ■ V0 d(f) t 1 d<p t 
q{T) = = aW P = ~2vr # • (A ' 6) 

From (A. 6) the representation (A.l) is rewritten in the Clebsh representation: 

B = V^xV(^-4 (A.7) 

The general form of the vector potential A is deduced from (A.7). Indeed, B = 
Vip x V0 + x = V x (-0V0 + qdVip) and the vector potential A is written in 
the Clebsh parametrization [44]: 

A = + q^)'§Vr + Vg, (A.8) 
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where g is a gauge function*. When the gauge d^g = is chosen, then the toroidal 
component = A ■ e$ is 

A^ = il>/R, (A.9) 

because the unit toroidal vector is = RV(j>. 

In section 2.1 another representation of the magnetic field is used. The B toroidal 
component is written as: 

^-V(p t xVtf = FV0, (A.10) 

where d^F = for axisymmetric systems and d$F = OJj when the condition Vr • J = 
is imposed on the plasma density current J = V x B/(4tt). Indeed, 

4vrVr • J = V • B x Vr = VF • V0 x Vr = d^F/^gJi^ = 0, (A.ll) 
because of the following equivalences: 

V • (VrZ> x V0) x Vr = V • (^'|Vr| 2 V0) = V(^'|Vr| 2 ) • V0 = 0. (A.12) 
(A.l) is now rewritten as 

B = V^/j x V0 + F(il>) V0. (A.13) 
The flux surface average of some function f(x) is defined as 

(/) (r) = / d 3 5;f(i)5(r - f) I [ d 3 S5(r - f). (A.14) 



In flux coordinates, (A. 5) allows to write (A.14) as follows: 

(/) (0 = / fy/g^tiWl f Vg^dW = f f f^l f f^. (A. 15) 

Appendix A. 2. Shafranov coordinates 

Shafranov coordinates are useful for axisymmetric toroidal equilibrium geometry 
characterized by nested flux surfaces with circular cross sections. The difference with 
the standard model [45] is on a relative shift of the centers of the circles corresponding 
to different flux surfaces: the Shafranov shift A(r). The map between the cylindrical 
coordinates (R, 0, Z) and the Shafranov coordinates (r, 9, <f>) is the following: 

' R = Rq- A(r) +rcos6» 

Z = rsin# (A. 16) 

$ = <fi, 

where R is the major radius of the magnetic axis and A(r) is normalized to give 
A(0) = 0. The center R g of the outermost magnetic flux surface when r = a is obtained 
if A(a) =Rq- R g . 

* It is worth noticing that in this representation A is a multivalued function. This is not a problem 
because it is consistent with the gauge invariance of B = V x A, because qip' (■d+2n)Wr = qip'tiVr — Vifit- 
jj This result can be obtained for a general poloidal angle 9 different from following the same steps 
of the sketched demonstration. 
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Figure Al. Displaced circular magnetic flux surfaces of Shafranov geometry. 



This geometry is qualitatively depicted in Figure Al where the directions of the unit 
orthogonal vectors e r , e$ and respect to the unit vectors e# = VR and ez = VZ, are 
shown. Once the V operator is applied to (A. 16), the following relations are obtained: 

e r = e» cos 9 + ez sin 9 = (1 — A' cos 9) Vr 

(A.17) 

cq = —eR sin 9 + ez cos 9 = A sin #Vr + rV#. 
Equations (A.17) can now be reversed requiring A' < 1: 

Vr = (1 - A'cos^rV 

, , 1 ( A -!8) 

rW = +e e - A sin ^(1 - A cos 9) e r . 

From (A. 18) and using the orthogonality property of the left handed basis e r ,ee,e0, it 
is computed the Jacobian: 

v ^7= (Vr x V0- V^ 1 = (1 - A'cos9)rR. (A.19) 

Once obtained the Jacobian, the expression for computing the flux surface average 
(A. 14) of a generic function / = f(r,9) is straightforward obtained: 

(/> (r) = / f^dB/f^dB = I j f ^Z^tn d9 ' < A - 2 °» 
because of the following integration: 

1 - A' cos 9)rR d9 = 2tt \(R - A)r - r 2 A'/ 2] . (A.21) 
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From (A. 13), the magnetic field is expressed as: 

B = Rd-L^f + F R e *- < A - 22 » 

The correspondence with the flux representation is obtained computing the flux poloidal 
angle d from the Shafranov poloidal angle 9. The relation $ = ■$(#,r) is obtained from 
the integration of old = devoid + d r ddr on a constant flux, or constant r, path. 

Equations (A. 13) and (A. 19) give 
_ . 5-W 5-W , B-V(f) , rF(l-A'cos0) , AOO , 

ded = Wve = — ^v^7 = -^r-V^I = ^ ( A - 2 3) 

and 

W = - 'I » \ /•" ' A '^° dO. (A.24) 
q^(R -A)J 1 + rcos9/(R - A) 1 ; 

From (A.24) and the 2-7T periodicity of the poloidal angles, ip' can be expressed as 

*' = -^A) X ^ (A - 25) 

where 

o-/ x 1 f 1 - A' cos # - 

X r = — <b = dO = 

2nl l + rcos9/(R -A) 

= 1 | 1 W r V fc ~V r | a A (2fc-l)H 
^27T^ 1 U -Ay' Uo-A^ 7 2A;!! ' 

When the inverse aspect ratio e = a/R g is little enough and A' = 0(e), the relations 
(A. 20), (A.24) and (A. 25) are truncated and respectively approximated by: 

(/) (r) = ^ / d0/{l + [r/(i? - A) - A] costf + 0(e 2 )], (A.27) 

= 9 - [r/{R - A) + A'] sin + £>(e 2 ), (A.28) 

*'-^A) + ^- < A - 29 » 

where it is assumed ■&{9 — 0, r) = 0. 

Moreover, the magnetic field is written as: 

rF[l + C(e 2 )l F lk , 

from which it follows the magnitude used in (26): 
F 



^2 



1 + 2?(i^ + ^ 



(A.31) 



Sometimes it is preferred the usage of geometrical coordinates when the poloidal 
sections of the flux surfaces are circular as the examined case. The transformation map 
can be obtained by (A. 16) with the following relations: 



R = R g + r g cos 6 g 



r g sm9 g (A.32) 



From the orbit theory to a GC parametric equilibrium distribution function 34 



The geometrical coordinates (r g , 9 g ) can be expressed in terms of the Shafranov 
coordinates (r, 9) by the following reletions: 

'r 2 = r 2 + 2rA Q (r) cos 9 + A 2 (r) 
g 9\ ) g \ ) A _ 33 

cot9 g = cot 9 + A g (r)/(r sin 9), 

with A g (r) = R — R g — A(r). When A g (r)/r <C 1 it is possible to approximate 
r ~ r g — A g (r g ) cos9 g . In this case it is straightforward to re-write the relations (A. 27, 
A. 29, A. 28, A. 30, A. 31) as function of (r g , 9 g ). As for example the poloidal flux expressed 
in geometrical coordinates ip(r) = ip g {r g ,9 g ) becomes: 



ip{r) = if; dr ~ / if; dr ~ 

/ r 3 
i>' dr - A g cos 9 g ip'(r g ) = ifj(r g ) - A g ip'(r g ) cos 9 g , 



(A.34) 



so that the Fourier representation of the equilibrium flux potential ip g ~ ip{j' g ) — 
A g ip'(r g ) cos9 g in geometrical coordinates involves only the first harmonic. It is worth 
noticing that (A.34) depends on the geometry, not on equilibrium constrains as can be 
the application of the Grad-Shafranov equation. Moreover, (A.34) fails on describing 
a system near to the magnetic axis where A g {r)/r ^ 1 does not occur because 
A g (0) = Ro — R g 7^ is assumed. 

The reader interested in the generalization of the Shafranov coordinates for 
describing axisymmetric plasmas subject to the Grad-Shafranov equation in shaped 
poloidal magnetic flux surface sections, can consult for example [46]. 
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